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Abstract. In this paper we study the apphcation of the Sobolev gradients technique to the prob- 
lem of minimizing several Schrodinger functionals related to timely and difficult nonlinear problems 
in Quantum Mechanics and Nonlinear Optics. We show that these gradients act as preconditioners 
over traditional choices of descent directions in minimization methods and show a computationally 
inexpensive way to obtain them using a discrete Fourier basis and a Fast Fourier Transform. We 
show that the Sobolev preconditioning provides a great convergence improvement over traditional 
techniques for finding solutions with minimal energy as well as stationary states and suggest a gen- 
eralization of the method using arbitrary linear operators. 
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1. Introduction. The observation of Nature reveals that many unforced con- 
tinuous systems tend to accommodate into stationary configurations, in which the 
distributions of mass, charge, velocity, etc, do not change throughout time. In the 
language of mathematical modeling all configurations are represented by points of a 
certain space of functions, ip{x.) G W, while the tendency of the system to lie into any 
of these states is given by a functional, E^ip) : W M., the energy, whose minima 
are precisely those stationary "states". For this reason it is possible to see many 
physical problems written as variational principles of the type "find ip € W such that 
E{ip) : W ^ M. achieves a minimum on W" . In most situations the functional to be 
minimized has a dependence on -0 of the form 



However, a complete analytical description of the minima of the functional E(ip) is 
usually not possible. In this paper wc will introduce several techniques for performing 
this study numerically, focusing on the minimization of subject to physical 

constraints. 

From a practical point of view, this problem is similar to that of finding the 
minima of a real function defined over a finite-dimensional space, such as R". First, a 
definition of derivative of the functional, Vi?(V'), must be chosen. If the domain of the 
functional W is equipped with some scalar product we may use the Frechet derivative 
which is given by a first order expansion of the functional around a function ip 



The critical points, tpc, are defined as the points where the first order variation 
of the functional vanishes for any perturbation 5. That is, the derivative vanishes in 
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a weak sense 

(1.3) (<5,v£;(^e)) = 0, ys. 

Just like in the finite-dimensional case it is possible to show that any minimum of 
the functional must be also a critical point. Thus a common approach is to solve Eq. 



(1.3) and verify a posteriori which solutions are actually minima of the fmictional. If 
W = L^(R") this procedure gives us the well known Euler-Lagrange equations of the 
problem, which is a partial derivatives equation (PDE) 

(1.4) |£_v.#L = o. 

dip oVip 

However, we have no guara ntee to reduce the complexity of the problem, as it is 
by no means trivial to solve Eq. (|l.4| ). Also we are likely to obtain more solutions than 
we actually need, since not only minima, but also maxima and saddle node points 
will satisfy the Lagrange equations. 

To avoid these problems some other methods are used which aim at finding the 
minima of the functional directly, constructing minimizing sequences, {ipi}, whose 
limit is a minimum of the functional: ■0c = linii^oo ipi- These methods will be dis- 
cussed in the following sections. 

The outline of this paper is as follows. In we recall the definition of Sobolev 
gradients as given in ^, ^. We derive a formal solution to the problem of finding 
these gradients which is based on the inversion of a positive hermitian operator. In §^ 
we derive an explicit expression for the Sobolev gradients in the trigonometric Fourier 
basis and comment on its implementation using Fast Fourier Transforms (FFT). 

In §11 and §|| we apply the previous tools to two physical problems. Using descent 
techniques with Sobolev gradients over Fourier spaces we will find the ground states 
of a Bose-Einstein condensate in a rotating magnetic trap and the excited states for 
coupled laser beams propagating through a nonlinear medium. Both physical systems 
are modeled by nonlinear equations of Schrodinger type and present difficulties when 
traditional minimization techniques are used. We comment on the great improvements 
that are achieved using Sobolev gradients. Finally in §^we summarize our results and 
offer some conclusions. 

2. Sobolev gradients. 

2.1. Direct solutions of the variational problem. There are two traditional 
approaches to the problem of finding the minima of a functional using a discrete basis. 
The first one expands the unknown solution using a Fourier basis {4>k\, ip = X]/c '^k4'k, 
and defines a new functional over the finite-dimensional space 

(2.1) E{{ck})^E(j2ck^ 

The functional is then minimized using methods which are well known from the the 
domain of finite-dimensional problems, e.g. Newton's method or nonlinear conjugate 
gradient. 

This procedure is quite straightforward and there is a huge amount of literature 



and tools which can be immediately applied to (2T). However, for some types of 
problems one has to deal with highly nonlinear algebraic equations with many of 
terms on each equation, something which is computationally too expensive to work 
with. 
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The second approach involves what is known as descent techniques. The idea 



is to manipulate the original functional (1.1) building an analytic equation for a 



minimizing trajectory in the target space W . This equation is then discretized and 
solved on a suitable basis. In the continuous steepest descent version, the trajectory 
tp{t) : M. ^ W is continuous and defined by a PDE which involves the gradient of the 
functional ( |l.2| ) 

(2.2) f = -VE. 

The discrete steepest descent technique is computationally cheaper since instead of 
requiring an integrator for the PDE it constructs a discrete succession of estimates to 
the minimum, {ipk+i = V'fc + AfeVi?('0fe)}, by locally minimizing E^tpk + AV-B) with 
respect to the real parameter A. 

In this paper we will only deal with descent techniques. The first and most 
important reason is that we will work with the definition of V E(ipk) trying to improve 
its convergence. As it was already shown in , this work pays off: a good choice of the 
gradient improves convergence by several orders of magnitude. The second motivation 
is that by focusing on the gradient our algorithms will be essentially independent 
on the descent method, which leaves space for further improvement. For instance 
one might apply these techniques to a nonlinear conjugate gradient method — which 
dynamically adjusts the search direction, dk = VE{^k), with an estimate that takes 
into account the history of the evolution — . Finally, by working with Eq. ( |2.2|) or its 
discrete version we avoid the complex nonlinearities that arise in other methods, and 
we will have a more ample choice of Fourier basis to work with. 

2.2. Ordinary gradients. It is customary in the literature to work in spaces 
which are equipped with an scalar product and its corresponding norm 



(2.3) = J 

(2.4) Mi^ ^ I ivr 



If one does so and works with Eq. (E2) then the formal definition of the gradient is 
Lagrange's one 

(2.5) Vi?(V') = ^-V- 



9V d (V7/i) ■ 

We will refer to this definition as the "ordinary" gradient, to distinguish it from the 
different definitions that we will derive below. 

2.3. Sobolev gradients. Following the ideas from j|] we will move our problem 
to a different space which is the Sobolev space of functions, such that ^ and its 
derivatives, Vi/", have a well defined L^-norm: 

(2.6) = {V'/V',VV' e i^}. 

This Sobolev space will also be equiped with a scalar product and a norm 

(2.7) (V, 4>)= I [V;(x)0(x) + V^(x) • V0(x)] 



(2.8) IIV'II'^ / [I^P + lVV-Plrf 
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To obtain a new explicit expression for the gradient of the functional in the 
Sobolev space we will follow a less rigorous derivation than in j^] . Performing a first 
order expansion of our functional around a trial state ip we obtain 



(2.9) Eii/j + e6) = E{iIj) + e {6,V6) { || ] + c.c. + 0{e^). 




We have to turn this expression into something like Eq. ( |1.2|) . This means that 
we have to find some 6 such that 



If we integrate by parts and impose that this equality be satisfied for all perturbations, 
i5, the problem has a formal solution which is given by a Lagrange equation 

.X , dE „ dE 
(2.11) (i_A)0=-y 



In consequence, our formal expression for the Sobolev gradient of E{ip) finally reads 
(2.12) VsE={l-l\y^VE. 

Here V sE stands for the Sobolev gradient, VE is the ordinary one, and (1 — A)^^ 
represents the inverse of a linear and strictly positive definite operator. 

3. Sobolev gradients on Discrete Fourier spaces. In the rest of the pa- 
per we will work with functions which are defined over a d-dimensional rectangu- 
lar volume = {x G ni[ai,6i]} with side lengths given by Li = hi — ai. We 
customarily define an orthogonal set of basis functions, 0„ = e*'^"^, over where 
{fc„ = 2^(fi,...,f^),n,ez} 

It is well-known that it is possible to expand any continuous function f{x) with 
periodic boundary conditions using this basis 

(3.1) /(x)= ^ /„0„(x) 

n— — c« 

where 

(3.2) /" = / Mx)fi^)- 

In the previous formula V ~ Hi Li is the volume of and arises because of the 
lack of normalization of the basis functions, a common practice which saves some 
computation time. 

To discretized the problem we will work within the set of functions sampled over 
a set of evenly spaced points from Q, {xn — {n\h\^ . . . , ridhd), rii = 0, . . . , A'^; — 1}. 
Here hi represents the spacing along the i-th dimension, n is a vector of non-negative 
integers and we denote a sampled function by an index, as in /„ = /(a;„). 

Due to this discretization our previous Fourier basis is now redundant. It can 
be reduced to a finite subset of functions which may represent any sampled function. 
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These functions are given by |fc„ = 2n (^j^, • ■ • , , rii = —Mi + 1, . . . , Mj| , where 

Mi — [^^^72] is the integer quotient of Ni divided by two. In this basis a sampled 
function is given by an expansion which reads 



(3.3) fm = y] fn4'niXm) 



and which makes use of the same coefficients as Eq. (3.2). 

The advantage of the finite Fourier basis over other approaches is that it provides 
an approximant of any function whose error is of order 0{Li/ N^Y where p is the max- 
imum differentiabifity of the sampled function. Furthermore, there is a numerically 
efficient method known as Fast Fourier Transform which allows one to compute the 
Fourier coefficients up from the sampled function, /„ — > f„i, and vice versa. 



Solving Eq. ( 2.12 ) numerically in a discrete Fourier basis is simple. Let us say 
that we have computed the ordinary gradient and that its sampled version has some 
Fourier coefficients 



(3.4) Vi;(x„) = e„i0„i(x„) 



and let us assume that there exists a certain solution to Eq. (2.12) and that it has 
another discrete Fourier expansion 

(3.5) Vs-B(x„) = ^ s,„0,„(a:„). 



Then by virtue of Eq. ( |2.12| ) 
(3.6) 



That is, in the sampled space the Sobolev gradient represents a preconditioning of the 
ordinary gradient, such that the most oscillating modes are more attenuated. Fur- 
thermore, due to this very simple expression computing the Sobolev preconditioning 
is computationally cheap and involves only minor changes to existing computer codes 
based on Fourier transforms. 

4. Applications to Quantum Mechanics. 

4.1. The problem. In this section we will apply the Sobolev gradients technique 
to a timely problem from Quantum Physics. The system that we will study is a dilute 
gas of bosonic atoms which are cooled down to ultralow temperatures at which their 
dynamics become synchronized. When the temperature is low enough, the gas or 
"condensate" may be described using a single complex wave function, ^{x,t) which 
is ruled by the so called Gross-Pitaevskii equation, a type of Nonlinear Schrodinger 
equation that for the Bose gas in a rotating trap with adimensional units reads 



(4.1) i5tV'(x,t) = 



-^A + V{^)+gm^,t)\^-nL, 



V'(x,t). 



Here g E M+,il E W, — i {xid-z — X2di) is an hermitian operator whose expected 
value (Lz) = J ^pLzip-, represents the angular momentum of the condensate along an 
axis of the trap. 
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There is conserved quantity associated to equation (4.1) which is called the energy 
functional of the condensate 



(4.2) 



F(x) + iglVp - nL, 



Our objective in this part of the work will be to find the solutions which are the 
minima of the energy subject to a restriction of the L^-norm 



(4.3) 



N. 



The particular value of N is imposed by the experimental conditions. Furthermore, 
without this restriction the absolute minimum of the energy is always reached at the 
trivial solution ijj = 0. 

In Quantum Mechanics the variational formulation of the problem is traditionally 
converted into a Lagrange equation (1.4), which is nothing but the Gross-Pitaveskii 
equation for so called stationary states. In short the word "stationary" refers to 
solutions of the type 



(4.4) 



V'(x, t) 



The pair {/i, (/)p(x)} satisfies a nonlinear eigenvalue problem 



(4.5) 



-iA + y(x)+.g|0^(x)|2-r!L, 



Due to the difficulty of solving directly the problem (4.5) we will try to find a direct 
solution to the variational problem. 

4.2. N um erical methods: imaginary time evolution. We can search the 
minima of (4.2) using descent techniques modified to account for the restriction on 
the norm (4.3). The first way to do this is to use a version of the continuous steepest 
descent which is known as imaginary time evolution [?]. We can summarize this 
method with the following set of equations 



(4.6) 

(4.7) 
(4.8) 



cr(x, r). 




Here we see that i/(x, r) evolves continuously maintaining a fixed L^-norm N and 
following the direction of decreasing energy given by Vi?. Indeed it is easy to show 
that £^ [E{i'{x,t))] < 0. Hence, the limit given by 

(4.9) (/)(x) = lim i^(x, t) 

r — >oo 

is at least a critical point of the energy, if not a minimum^ 



^As it is common with these local minimization procedures, it remains the problem that iterations 
may be trapped in a critical point which is not a minimum. Linear stability analysis may then be 
applied to check the validity of the solution. 
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From a practical point of view, the simplest way to find the minimizer using 
imaginary time evo lution is to repeatedly integrate Eq. (4.7) for a very short time, 
At, apply Eq. (4.6) for the newly found (t(x, t + At) and use the new estimate for v to 
redefine a = v and repeat the procedure until convergence. This way one avoids the 
problem that according to Eq.(4.7) the norm of a may grow indefinitely. The same 
consideration applies to the remaining methods that we will present here. 

Although the method from Eq. (4.7) was derived using an ordinary gradient, 
nothing prevents us from applying our Sobolev preconditioning and all results should 
be still valid. If we do so, our new equations are 



(4.10) 

(4.11) 
(4.12) 



cr(x,T), 




y(x)+5|i.|2-r!L, 



and the critical point is still given by Eq. 

4.3. Numerical methods: minimization of the free energy. While the 
imaginary time evolution is easy to understand and to implement, the fact that it 
performs the descent over a path of functions with a certain norm makes it too re- 
stricted and sometimes too slow. A different approach is to define a new functional 
called the /ree energy with a Lagrange multiplier that takes care of the restriction on 
the norm. In Quantum Mechanics this free energy is usually defined as 



(4.13) 



because it preserves the linear form of the equations. Since our equations are already 
nonlinear we define a free energy functional more conveniently as 



(4.14) 



F{^) = E{^) + i (iV(^) - A)' 



First and most important it is not difficult to show that any absolute or relative 
minimum of -F'(?A) is also a minimum of ^^("0) subject to Eq. ( |4.3| ) with a nonlinear 
eigenvalue given by /i — N{ip) — A. 

Secondly, as we will prove in Appendix F{'ip) must have at least one finite 
norm minimum, something which cannot be easily assured for Fqm- 

The practical advantage of our new functional consists in that fixing A and fl we 
can perform a continuous descent over the whole domain of F^ip) without renormaliz- 
ing the solution on each iteration — i.e., the search space is larger. The new equation 
that we must integrate is thus 



(4.15) 



97 ("'"^ 



-VF{iy) ^ - 



--A + y(x)+g|z/|2-r!L, 



Using a Sobolev preconditioning Eq. (4.15) becomes 



(4.16) ^(x,t) = -VsF(z.) = -(1-A) 



In both cases the actual solution is still given by the limit of Eq. 
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4.4. Numerical results. Up to this point we have shown four different numer- 



ical methods, two of them incorporating a Sobolev preconditioning [Eqs. (4.11) and 
( 4.16 )] and two without it [Eqs. (4/7) and ( 4.15| )]. We have compared the efficiencies 



of these methods for several test situations. The details of our study are as follows. 

All methods have been implemented using the discrete Fourier basis mentioned 
above. This applies both to the calculation of derivatives and to the application of the 
Sobolev preconditioning. We restricted our simulations to two-dimensional problems 
with radially symmetric potential, V{x) — jjxp, over a grid of 128x128 points. 
Practical experience with more complex problems shows that our results generalize 
to higher dimensionality and denser grids. 

Due to the requirements of the method, both variants of imaginary time evolu- 



tion [Eq. (4.7) and ( 4.11 )] are integrated using a Runge-Kutta-Fehlberg, where the 
tolerance is adapted as the solution converges to its target. On the other hand, in- 
stead of performing a continuous descent for both variants of the free energy descent 



[Eq. (4.15) and (4.16)], we found it more convenient and faster to perform a discrete 
steepest descent. 

All programs have been implemented using the tensor- algebra environment called 
Yorick ]|6| , an interpreted environment which is capable of fast numerical computations 
and which can be equipped with the FFTW library [Q . Execution times are given for 
a Digital Personal workstation 500au, but the programs run equally well on modest 
personal computers with Pentium-II processors and less than 64 Mb of memory. 

As test case we have considered three situations. Case A is the simplest one 
and corresponds to a stationary trap 51 = 0, an intense nonlinearity g = 100 and 
starting the minimization with a radially symmetric Gaussian of unit width, that is, 
"i/iQ (x. e"'"' , a shape which is similar to that of the true ground state. 

Cases B and C involve a rotating condensate with fl ~ 0.6 and g = 100. Under 
these conditions the functional achieves both an absolute minimum, which is the same 
that found in case A, and a local minimum. The local minimum is a solution of vortex 
type, a topological defect, whose behavior near zero is tp (x {xi + ix2)/\x.\'^- 

For this reason we designed two test cases with different initial conditions. Case B 
starts with a Gaussian profile with a centered vortex or ipo (x [xle^l'^l (xi -|-ia;2)/|x|^. 
In this case all methods are trapped on the local minimum with the centered vortex. 

Finally, a third set of simulations, case C, uses the same parameters {51 — 0.6, g = 
100} but starts from a nonsymmetric initial state, ipo oc |x|e~l''l ((xi — yi) + i{x2 — 
j/i))/|x — y|^ which is close to the vortex solution but belongs to the basin of attraction 
of the ground state. Here all minimization methods require more computational work 
since they must find a path out of the local minimum, and it is precisely in this case 
where the differences between methods are best shown. 



In Table 4.1 we summarize the results of the simulations. In case A it is apparent 
that the Sobolev preconditioning has a positive influence over convergence, with an 
astonishing result of 55 steps for the steepest descent with free energy. A similar 
behavior is found in case B. 

In case C, the Sobolev preconditioning enhances convergence speed by two orders 
of magnitude. An intuitive explanation of why the steepest descent with a Sobolev 
gradient takes less steps to converge will be discussed in detail in Appendix 

5. Applications to Nonlinear Optics. 

5.1. The model. In this section we will consider a model for a pair of inco- 
herently interacting light beams. To be precise we will study the light field of each 
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Fig. 4.1. Evolution of error, £2 = \\ip — ^exact\\2, through different minimization processes for 

continuous steepest descent with Sobolev preconditioning (lower solid line) and without it (dashed 
line), and for imaginary time evolution with Sobolev preconditioning (upper solid line) and without 
it (dotted line). Plots (a,) to (c) correspond respectively to the cases A,B, and C described in the 
text. Both axes, error and number of iterations, are in a logarithmic scale. 
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ITS 


FE 


FES 


Case A 


Iterations 


1320 


945 


2850 


55 


Time (s) 


416 


371 


285 


13 


Case B 


Iterations 


1630 


615 


3210 


320 


Time (s) 


468 


242 


75 
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Case C 


Iterations 


64195 


2665 


108505 


1455 


Time (s) 


19863 


1165 


10861 


168 



Table 4.1 

Iterations and computation time for each minimization method: imaginary time without (IT) 
and with (ITS) Sobolev preconditioning and free energy without (FE) and with (FES) Sobolev pre- 
conditioning. Shown are results for the initial data described in the text (cases A, B and C). 



beam, u(x, and w{x,t), propagating through a weakly nonlinear saturable optical 
medium. This system may be modeled by the Cauchy problem 



(5.1) 
(5.2) 



idfU 
idtw • 



-Au + 



-Aw + 



1 + k(|«|2 + |«,|2)' 
w 

1 + k(|u|2 + |w|2)' 



for the complex functions u,w : M.^ x IR+ — > C, which vanish at infinity and satisfy 
the initial data u(x, 0) = uo{x.) and w(x, 0) = wo(x). Here n e R+, —A = ^3 + ^ 
is the Laplacian operator which accounts for the diffraction of light and the nonlinear 
term (l + jup + juip)"^ models the saturable interaction among the beams. 
Let us define the two component vector 



(5.3) 



[/(x,t) 



u(x, t) 
w(x, t) 



then, the energy functional for the system reads 



(5.4) 



E{U) = J [-;7+Af/ + G(|C/n 
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where G{p) — (hi (1 + up) — np). The analysis of this section may be generalized 
to more general nonlinearities with only minor changes to G{p). 

5.2. Stationary solutions. We are interested on stationary solutions, which 
are of the form 

(5.5) U{^,t) = "'q"* f^^, ) C/(x) = e'*"[/(x). 

The equations for the stationary solutions now are of elliptic type 

(5.6) MU = -AU + G'{\U\'^)U. 

with zero Dirichlet boundary conditions at infinity. Again this formulation poses a 
nonlinear eigenvalue problem for the pair {M, [/}. 

The stationary solutions are critical points of the energy functional subject to a 
constraint on the L^-norm of each component. That is, defining 

(5.7) ^" " / 

(5.8) N^^ I 



the first order variation of the energy around C/q for fixed and must be zero 

5E 



= 0. 



The particular fixed values of {Nu, Nyj} represent the total intensity of each beam of 
light. 

5.3. Ground states. In principle there are many different stationary solutions 



of Eq. (5.6). However, if we focus on the ground states (stationary solutions of 
minimal energy) we can apply some of the methods that we mentioned in with 
minor modifications to account for the higher dimensionality of the problem. For 
instance, one may define a free energy functional 

(5.10) F{U) = E{U) + \{N^- Kf + \{N^- , 

and minimize it using a discrete steepest descent with Sobolev preconditioning. 

By performing this minimization using different parameters {A„, X^} one obtains 



nodeless localized solutions for u and w as shown in Fig. 5.1. The precise values of 
Xu and Xw determine the norm of each component. 

Let us remark that, up to a global factor, the ground state has the same shape 
in both envelopes. That is, 

(5.11) uoi^)=N^pi\^\,N), 

wo{^)^n^T^p{\kIn), vee[o,i]. 

The common shape p depends on the total intensity N — + N^, and it is the one 
that fixes the values of and pw ■ For this reason if we had chosen the traditional 
definition of the free energy Fqm{U) ~ E{U) + puNu + PwN^ we would have found 
an infinite degeneracy which prevents convergence to the desired values of {N^ Nw\. 
This problem is removed by the use of our nonlinear Lagrange multipliers and the 
{Xu,Xiu\ parameters. 
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Fig. 5.1. (a) Density profile of a ground state for Nu = Nw = 30. (b) Norms Nu = Nw = N, 
as a function of the nonlinear Lagrange multipliers Xu = Xw = X. 



5.4. Excited states. In the field of guided light waves there has been a great 



interest on the properties of solutions of Eq. (5.6) which are not ground states, the so 
called excited states. Minimization methods based on the energy functional [See §^ 
cannot be applied to this task since excited states are not necessarily local minima of 
the energy. Instead, a variational principle somehow equivalent to Eq. ( |5.6D must be 
defined. 



Let us rewrite Eq. (5.6) as the application of a nonlinear operator 

(5.12) /(C/)ee [-A-A/ + /(|C/|2)][/ = 0, 
define the error functional 

(5.13) F{U) = I f{U)^f{U) > 

and take our variational principle to be "find Uo such that F reaches an absolute 
minimum F{Uo) = 0." 

With this principle and a Sobolev preconditioning our descent technique becomes 

(5.14) ^^igEll) ^_VsF(z.) = -(l-A)-iVi^, 

and it may be proven that Uq{x.) = liniT-^oo '^(x, r). The ordinary gradient V F reads 

(5.15) G^{-A + I{\U\^))U, 

(5.16) VF = (-A + I{\U\^))G + I'i\U\^) {U'^G + G'^U) U, 

Its calculation using a discrete Fourier basis is straightforward. 

There are several advantages of this approach. The first one is that F{U) makes 
no distinction between ground and excited states: any stationary state with the right 
eigenvalues fii — Ala is a local minimum of this new functional. The second one is 
that we do not need to add any Lagrange multipliers to F{U) since they arc already 
present in the M operators. Finally we expect that the number of minima of F{U) be 
discrete and separated so as to avoid problems with descent methods being trapped 
on critical points that are not minima. Indeed, it is very easy to know when this 
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Fig. 5.2. (a) Vortex-mode vector solitons and (b) their norms Nu, N^, Nt^tal = + N^, 
as a function of fiu for fim = 1. (d) Dipole-mode vector solitons and (c) their norms Nu, N^i, 
Ntotal = Nu + Nui, as a function of fiu for fi^i = 1. 
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Fig. 5.3. Different unstable multi-solitonic configurations arising from Eq. 
change of the initial conditions for fixed fiu = —1, Atui = —0.3, k = 0.5. 



(5.14) with a 



accidental trapping happens since any absolute minimum of F must also be a zero of 
it F{Uo) = 0. 

A remarkable feature of the method that we have outlined above is that to dis- 
tinguish among the different excited states we have to change both the ad-hoc eigen- 
values, {/Lt„,/iu)}, and the initial data of the minimization method. 

We have concentrated on two types of excited states of particular interest for 
applications: the first one is called vortex vector soliton its features being summa- 
rized in Fig. |3.2K a-b). Choosing different initial data for the minimization process we 
obtain an second type of excited states which are called dipole-mode vector solitons 



An example of these asymmetric solutions is depicted in Fig. 5.2(d). Depending 
on the parameters ^u,w we find different norms of the solution. The fact that our 
method allows the finding of these nonsymmetric stationary states is very interesting 
since conventional approaches to the problem have severe difficulties. In fact, the 
application of the Sobolev preconditioning not only greatly enhances the convergence 
but it is necessary to get convergence. 

Nevertheles s th e method works equally well for more complicated stationary so- 
lutions. 



In Fig. 



5.3 we show several exotic solutions that are obtained by solving Eq. 
( 5.14 ) with different initial conditions. All of those states are dynamically unstable 
and could have not been found with traditional minimization methods. 
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Fig. 5.4. Dependence of the error on the number of iterations when looking for a (a) 
vortex-mode vector soliton or a (b) dipole-mode vector soliton. We show results for a grid with 
32x32 points (solid), 64 x 64 points (dotted) and 64x64 points starting from an interpolation of the 
32x32 solution (dashed). 



5.5. Performance and grid refinement. In this subsection we want to an- 



alyze the performance of the method from §5.4 and to introduce a multigrid-Hke 
technique to improve convergence rates while looking for more accurate sohitions. 

The idea of the multigrid technique is to use solutions from coarse-grain grids to 
calculate better approximations on finer grids Roughly the algorithm consists 

in setting a coarse-grain initial data, solving the equation or the variational principle 
with this initial data until the error is small enough, interpolating the result over 
a finer grid and iterating using this new grid until both the error and the spatial 
discretization of the solution are the ones we desire. 

Since we are already working with Fourier modes over discrete grids the logi- 
cal choice for our algorithm is indeed Fourier interpolation. The idea is to use the 
expansion from Eq. ( p.3[) outside of the original grid, that is, 

(5.17) V^"""^ (x) = J2 Cne*"'', Vx e n. 

n 

This approximation is then discretized over a finer grid, resulting an expansion with 
a larger number of modes 

(5.18) Vm""^ = J2 c("^'")e'''""""^", 

n 

which may be used as the starting point for further iteration. 

In our case we have used only two different grids, one with 32x32 points {coarse 
grid) and a other with 64x64 points (the fine grid). We have used two different initial 
data for our minimization method: either a pair of Hermite modes which resemble 
the desired shape, or the solution of the coarse grid interpolated over the finer grid. 



As both the evolution of the error in Fig. 5.4 and the computation times in Table 



5.1 show, interpolation saves a significant amount of time. Indeed, the interpolated 



solution with only 32 x 32 Fourier modes is already a good approximation for the fine 



jrid, as it shows the small change of the error in Fig. 5.4 



6. Conclusions. In this paper we have shown a numerically efficient way to 
improve the convergence of several minimization methods using the so called Sobolev 
gradients and applied it to different problems which involve Nonlinear Schrodinger 
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Initial data type 


Herm] 
for vort 


te modes 
ex solitons. 


Interpolated 
solution 
for vortex 


Hermi 
for dipo 


te modes 
le solitons. 


Interpolated 
solution 
for dipole 


Grid 


32x32 


64x64 


64x64 


32x32 


64x64 


64x64 


Iterations 


540 


2101 


812 


496 


1206 


496 


Time (s) 


106 


2006 


757 


60 


1141 


466 



Table 5.1 

Results for the search of stationary solutions of Eq. (5.1)-(5.i), with and without grid refine- 
ment. Shown are the results for different initial data. The initial data named "interpolated solution" 
corresponds to taking as initial data the approximated solution on the 32 X 32 grid. 



Equations. We have also proven that these gradients represent a preconditioning over 
the traditional definition of gradients on . In Appendix B we suggest a generaliza- 
tion of this method to different vector spaces and scalar products. 

We have presented two different methods for solving our minimization problems: 
a traditional one, the imaginary time evolution^ and a new one, the minimization of 
a nonlinear free energy. Both methods have been shown to be suitable for a Sobolev 
preconditioning, giving us two new methods that we call preconditioned imaginary 
time evolution and a preconditioned free energy. 

We have implemented all four methods using a discrete Fourier basis and a Fast 
Fourier Transform. With these tools we have shown that the Sobolev preconditioning 
becomes an inexpensive additional step over existing methods. The four resulting 
solvers have been applied to several realistic problems and in all tests the nonlinear free 
energy with the Sobolev preconditioning showed the best convergence rates. Indeed 
the Sobolev preconditioning has an important effect on convergence which can be as 
good as gaining two orders of magnitude over the traditional techniques. Furthermore, 
opposite to what happens with finite differences ||^ the preconditioning may be applied 
without significant computational cost. 

We have derived a new method for finding excited states of coupled nonlinear 
Schrodinger equations. This method introduces a new variational principle which is 
not based on an energy functional but on finding the zeros of a nonlinear operator 
which corresponds to the equation to be solved. We have also shown how to improve 
convergence using Fourier interpolation and a two-grid method as a source for better 
initial approximations of the iterative method. This approach may be extended to 
more sophisticated multigrid methods. 

Appendix A. Existence of minimizers for the nonlinear free energy 
functional.. 

In this appendix we want to prove the existence of minimizers for the free energy 



(4.14). Let us write the free energy functional in the following form 



(A.l) F{^) = / Mn^d'^x + I f + - Xf 



2'^' 2 



where N — J |^p(i"r and 



(A.2) An = -^A + V{^)-nL,, G [0, 1) 
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is a positive hermitian operator Aq > ^min > 0. 

Let us also define the following spaces of functions defined over R" . We will work 
in LP spaces 

(A.3) LP = i^ip : ^ C / U\\p (^J |^/.(x)|Pd"a:^ < +oo| . 

Also of interest will be the space of functions over which An is well defined 
(A.4) Hn = \^:R^^C/ Mn := ( f ^(x)AoV(x)) ^ < +00] . 



Let us remark that the dimensionality of the space, is important, since for n < 2 
it is easy to show that Hn C L*. 

With the preceding notation we will state the following relevant theorem: 



Theorem A.l. The free energy functional F given by Eq. (A.l) has at least one 
absolute minimum in the set given by the inequalities 

o<\mi<\\mi<\{\ 



Proof. The first step of the proof will be to show that the domain of F is indeed 
the whole space Hq. Using the positivity of the An operator we show that 

(A.5) IIVII2 < \/aw||V'||o, 

which means that Hn C L^. We need Sobolev's inequality ||ll| 

(A.6) IIV^IU < IIV-ll^-'llvv^ll^ 

where d — n/2^n/k, n is the dimensionality of the space and for us d — 1/2. Applying 
this inequality we obtain the following bound 



(A.7) ii^iu < vuhWW2 < i^umn, 

which means that Hq C L^. Since || • II2, || ■ ||4, || • ||n < +00 inside we conclude 
that F as given by Eq. (A.l) is well defined over the whole space. 

F is also continuous. To prove it let us rewrite the free energy functional as 

(A.8) i.(V,) = ||^||2 +|||^||4+1(||^||2_^)2^ 

and using the bounds ( [A.5| ) and ( |A.7| ) it is easy to show that 

(A.9) \F{^P)-F{0\<6{e), HV' - ^||o < e 

where the constant (5(e) is given by e and HV'IU- 
We will also need to show that F is coercive 

(A.IO) lim > a > 0. 



16 



J. J. GARCIA-RIPOLL AND V. M. PEREZ-GARCIA 



Using Eq. ( |A.S| ) one may show that indeed 
(A.ll) ^ > IIV'll 



and thus the quotient tends to infinity as the norm grows. The continuity and co- 
ercitivity of F ahow to use theorem 1 (section 1.2) of |lj], which states the existence 
of at least one minimizer {inf F{u) : u G X} of F provided it is a weakly lower 
semicontinuous and coercive application F : X — > R. 

So we know that there is at least one minimum, but we do not know where to 
look for it. Let us now show how the A parameter allows us to select different targets 
for the minimization problem. To do so we define a real one-dimensional function for 
any given direction G i?o 

(A.12) f{k) = F{k^). 

This function is nothing but a polynomial over k 

(A.13) f{k) - kNa^, + k^N^"^ + \{kN~ \f , 

where — ^ ipAip/N and ~ f / I'/^l*/-^^ constants that depend only on the 
precise direction ip. 

By differentiating the polynomial and imposing k = we find that /'(O) < at 
the origin for all possible directions. This means that, as we mentioned above, our 
Lagrange "multiplier" A allows us to avoid the useless solution, -0 = 0. 

Furthermore, we can restrict the location of the minimizer to a surface of a certain 
norm. By differentiating /(fc) we reach 

(A.14) (a^ - A) + Nk{u.^ + 1) = 0. 

This equation has a single solution which gives us the norm of minimal energy along 
the ^/i-direction 



(A.15) NrmnW = max <j 0, }<X-^^ 



rmn 5 



+ 1 

a value which is bounded above by our Lagrange multiplier A 



From Eq. ( A.14 ) it also follows that for any absolute minimum of the functional, 
ipmim the expected value of the operator must be bounded by the norm 



(A.16) J i^minAtprnin = HV'mmllfl < XN{tpmin)- 

Otherwise the trivial solution ip = would have less energy than ipmin- Wc can thus, 
instead of working with an unknown surface, delimit the location of the minimum to 
a set which is given by two inequalities, ( A.15| ) and ( |A.16 ) 



(A.17) W^i^peHn : < A^(0) < A - < AiV(V')}. 

□ 

There are several practical consequences of this theorem. First, it states that the 
problem of finding the minima of E{ip) subject to fixed norm has one solution, i.e.. 
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there exists at least one ground state. Second, but equally important, it proves our 
Lagrange penalizer ^{N — /i)^ to be specially well suited for this problem since it 
avoids both the useless solution ip = and those with too large norm. And finally it 
gives us a bounded in which the minimizer must be. Indeed we can extend this result 
by proving that we can restrict our search space to a compact superset of W. 

Theorem A. 2. Any absolute minimum, ip, of the functional F given by Eq. 
^■i ) lays inside the compact set of 

U = {ijeL^ : IIVII^ < A||^||2 < A(A-l)}. 



Proof. For the type of spaces that we work with, a set is compact iff it is closed 
and we can build an e-net for any positive number e. The e-net is a finite set if^ = 
{vi, . . . such that for each ip £ U there is an element v £ verifying — < 
e. Thus compactness is equivalent to the possibility of building an arbitrarily good 
approximation of our minimizer using a finite but sufficiently large basis of functions. 

It is evident that the set U is closed in the subspace Hq of L^. Let {«„} G C7 be 
a convergent sequence and let u be their limit. Since for each element of the sequence 

(A.18) ||u„||n< VA||u„||2<A(A-1), 

it is also obvious that 

(A.19) ||w||a< VA||ii||2< A(A-l), 

also thus the limit belongs to U. 

The compactness of the closed set U essentially follows from the fact that the 
eigenstates of Aq form a complete basis of Hq, and that the eigenvalues of An form 
a monotonously growing unbounded set of positive numbers. These eigenstates are of 
the form 

(A.20) ^„,^^^(|x|)::l±^e-l-l^/^ 

where P^^ are Laguerre's generalized polynomials. And the corresponding eigenvalues 
are 

(A.21) An(j)n,i = fin,i(l)n.i = {2n+il-n)l + l)(j>n,i, n,/eNU{0}. 

Let us choose any natural number k such that fc + 1 > A. We can split the whole 
space as a direct sum ~ H^^^ ^fc+i where 

(A.22) = \m{<j)„^i : j <2n + l < k}. 

The important point is that since i?Q^^ is isomorph to R™ for some natural number 
m, and Uk = U (1 Hq~^^ lays inside a compact m-dimensional ball of radius — 1, 
then we can find any e-net for Uk- Furthermore, by separating 

(A.23) = + i^b, i^a e Uk,iJb e H^^„ 

and using the definition of Uk we show that the projection of -ip outside of Uk can be 
made arbitrarily small 

(A.24) WMl < j^N. 
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A direct consequence of this is that for k > XN/e, a e-net of [7^ is also an e-net of U, 
which proves the compactness. 

Finally, by inspecting the eigenvalues of A^i and using E_q. (|A.15| - |A.16D it is not 
dilBcult to see that the absolute minimum of F must lay in U. 
□ 

Appendix B. Extending the Sobolev gradients. 

We have shown that redefining the gradient turns out to be a kind of precondi- 
tioning over the original choice of the direction of descent. Let us assume that our 
functional has the following form 

(B.l) Ei,p) = J v;a^ + /(ivp,x), 

where A is a non-negative hermitian operator that may involve some derivatives. Let 
us also assume that H is a suitable space equipped with the following scalar product 

(B.2) (V^,</)) = U{l + A)cj,, 



which is indeed a scalar product because A\^ > 0. 

Let us rewrite the energy functional in the following way 

(B.3) = (V-, V-) + J /(I^P, x) - l^p. 

The Sobolev gradient in H is 

(B.4) WAE^i^ + {l + A)-'[dif~i^], 

while the so called ordinary gradient is WE = Aip + dif. Hence the preconditioning 
nature of the method is recovered 

(B.5) VaE = {1 + A)~^VE. 

We can thus think that the new choice of the scalar product aims at making the 
linear part of the energy functional close to some quadratic form, [ij^^Bip), such that 
the new operator B is almost the unity. In our example, indeed, = tp. We that 
this preconditioning will both enhance the directions of decreasing energy and also 
have a smoothing effect on the nonlinear part. However, the problem is complicated 
and we know of no proof that these arguments be of general applicability. 



REFERENCES 

[1] .J. M. Sanz-Serna, Fourier techniques in numerical methods for evolutionary problems, in 
"Third Granada Lectures on Computational Physics^', Proc. of the III Granada Seminar 
on Computational Physics of 1994, P. L. Garrido, J. Marro, Springer- Verlag, Berlin, 1995, 
pp. 145-200. 

[2] ,J. W. Neuberger, R. J. Renka, Sobolev Gradients and the Ginzburg- Landau functional, SIAM 

J. Sci. Comput., 20(1998) pp. 582-590. 
[3] J. W. Neuberger, Sobolev Gradients and Differential Equations, Springer- Verlag, Berlin, 1997. 
[4] .J. J. Garcia- RiPOLL, V. M. Perez-Garci'a, Stability of vortices in inhomogeneous Bose- 

Einstein condensates subject to rotation: a three dimensional analysis, Phys. Rev. A, 

60(1999), pp. 4864-4874. 

[5] J. NOCEDAL, S. ,]. Wright, "Numerical Optimization" , Springer Series in Operations Research, 
Springer, New York, 1999. 



OPTIMIZING SCHRODINGER FUNCTIONALS WITH SOBOLEV GRADIENTS 



19 



[6] The Yorick environment for scientific comp uting is available at the electronic address ftp: //ftp- 



icf . llnl.gov /pub /Yorick / yorick-ad. html . 

[7J JVlATTEO 1^'RIGO, STEVEN G. JOHNSON, FFTW: An Adaptive Software Architecture for the FFT, 
in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Pro- 
cessing, IEEE, Seattle, WA, 1998, pp. 1381-1384. See also 

[8] Z. H. MusSLiMANi, M. Segev, D. N. Christodoulides, M. Soljacic, Composite Multihump 
Vector Solitons Carrying Topological Charge Phys. Rev. Lett., 84(1999), pp. 1164—1167. 

[9] ,1. J. Garci'a-Ripoll, V. M. Perez-Garci'a, E. A. Ostrovskaya, Y. S. Kivshar, Dipole-mode 
vector solitons, Phys. Rev. Lett., 85(2000), pp. 82-85. 

[10] W. L. Briggs, V. E. Henson, S. F. McCormick, "A multigrid tutoriat' , SIAM (2000). 

[11] L. Hormander, "The analysis of linear partial differential operators i". Springer, Berlin (1983). 

[12] B. DacorOGNA, "Direct Methods in the Calculus of Variations" , Springer- Verlag, Berlin (1989). 



